# This file estimates longterm effects of twinning for fathers who are 20 years older 

# install.packages("dplyr")
# install.packages("date")
# install.packages("AER")

library(dplyr)
library(date)
library(AER)

# Load data
load("fathers_old_20.rdata")

# Run OLS models for any children 

ols_dat1 <- 
  fathers %>% 
  mutate(children = ! is.na(age_oldest)) %>%
  filter(! is.na(stemt)) %>%
  mutate(vote = stemt == "Stemte")

ols_dat1 <-
  within(ols_dat1, {
    age        <- ntile(date, 25)
  })

ols_mod1 <- lm(vote ~ children + as.factor(age), data = ols_dat1)

# Filter of if age of oldest is unkwown (select parents )
fathers <-
  fathers %>%
  filter(!is.na(age_oldest)) 

# split data file to each election 
# 2009 election
# Subset and create vote variable
# filter of parents with children born after the election

fathers09 <- 
  fathers %>%
  filter(! is.na(stemte_2009) & first09 > 0) %>%
  filter(age_oldest <= (14565 - 20*365 - 5)) %>%
  mutate(vote = stemte_2009 == "Stemte",
         no_children = no_children_09)

#create age dummies for fathers 
#dichotomize twin first
fathers09 <-
  within(fathers09, {
    age        <- ntile(date, 25)
    twin_first <- twin_first > 0
  })

# create dummies for age of children grouped by age of father dummies
fathers09 <-
  fathers09 %>%
  group_by(age) %>%
  mutate(child_age = ntile(age_oldest, 4))

## Repeat for 13
fathers13 <- 
  fathers %>%
  filter(! is.na(stemt) & first13 > 0) %>%
  filter(age_oldest <= 16028 - 20*365 - 5) %>%
  mutate(vote = stemt == "Stemte",
         no_children = no_children_13)

fathers13 <-
  within(fathers13, {
    age        <- ntile(date, 25)
    twin_first <- twin_first > 0
  })

fathers13 <-
  fathers13 %>%
  group_by(age) %>%
  mutate(child_age = ntile(age_oldest, 4))

## Repeat for 14
fathers14 <- 
  fathers %>%
  filter(! is.na(ep_stemt) & first13 > 0) %>%
  filter(age_oldest <= 16215 - 20*365 - 5) %>%
  mutate(vote = ep_stemt,
         no_children = no_children_14)

fathers14 <-
  within(fathers14, {
    age        <- ntile(date, 25)
    twin_first <- twin_first > 0
  })

fathers14 <-
  fathers14 %>%
  group_by(age) %>%
  mutate(child_age = ntile(age_oldest, 4))

# Run models 
# Simple models. Biased but used for turnout of "control" group
mod09_r <- lm(vote ~ twin_first, data = fathers09)
mod13_r <- lm(vote ~ twin_first, data = fathers13)
mod14_r <- lm(vote ~ twin_first, data = fathers14)

# Models with fixed effect for fathers and children's age
mod09 <- lm(vote ~ twin_first + as.factor(age) * as.factor(child_age), data = fathers09)
mod13 <- lm(vote ~ twin_first + as.factor(age) * as.factor(child_age), data = fathers13)
mod14 <- lm(vote ~ twin_first + as.factor(age) * as.factor(child_age), data = fathers14)


## Instrument number of children with first child being twin
mod09_iv  <- ivreg(vote ~ no_children + as.factor(age) * as.factor(child_age) | 
                     twin_first + as.factor(age) * as.factor(child_age), data = fathers09)
mod13_iv  <- ivreg(vote ~ no_children + as.factor(age) * as.factor(child_age) | 
                     twin_first + as.factor(age) * as.factor(child_age), data = fathers13)
mod14_iv  <- ivreg(vote ~ no_children + as.factor(age) * as.factor(child_age) | 
                     twin_first + as.factor(age) * as.factor(child_age), data = fathers14)

## Run first-stage models 
mod09_fs <- lm(no_children ~ twin_first + as.factor(age) * as.factor(child_age), data = fathers09)
mod13_fs <- lm(no_children ~ twin_first + as.factor(age) * as.factor(child_age), data = fathers13)
mod14_fs <- lm(no_children ~ twin_first + as.factor(age) * as.factor(child_age), data = fathers14)

## Run simple OLS for relationship between no of children and turnout controlling for age fixed effects 

fathers13 <-
  within(fathers13, {
    no_children_ind <- no_children
    no_children_ind[no_children_ind > 5] <- 5
  })

ols_mod2 <- 
  lm(vote ~ as.factor(no_children_ind) + as.factor(age) * as.factor(child_age), data = fathers13)

# Compile results 

table_ols_13 <- 
  rbind(summary(ols_mod1)$coefficients[1:2,],
        summary(ols_mod2)$coefficients[2:5,])

table_dad_09 <-
  rbind(summary(mod09_r)$coefficients,
        summary(mod09)$coefficients[2,],
        summary(mod09_iv)$coefficients[2,])

table_dad_13 <-
  rbind(summary(mod13_r)$coefficients,
        summary(mod13)$coefficients[2,],
        summary(mod13_iv)$coefficients[2,])

table_dad_14 <-
  rbind(summary(mod14_r)$coefficients,
        summary(mod14)$coefficients[2,],
        summary(mod14_iv)$coefficients[2,])

ns <- rbind(
  with(fathers09, table(twin_first)),
  with(fathers13, table(twin_first)),
  with(fathers14, table(twin_first)),
  with(ols_dat1, table(children))
)

tab_first_stage <- 
  rbind(summary(mod09_fs)$coefficients[2,],
        summary(mod13_fs)$coefficients[2,],
        summary(mod14_fs)$coefficients[2,])

outtab_dad <- 
  list(table_ols_13,
       table_dad_09,
       table_dad_13,
       table_dad_14,
       ns,
       with(fathers13, table(no_children_ind)),
       tab_first_stage)


save(outtab_dad, file = "father_twin_old_20_res.rdata")

